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Abstract 

Using Heavy Quark Effective Theory with non-perturbatively determined parameters in 
a quenched lattice calculation, we evaluate the splittings between the ground state and 
the first two radially excited states of the Bg system at static order. We also determine 
the splitting between first excited and ground state, and between the B* and Bs ground 
states to order l/mb- The Generalized Eigenvalue Problem and the use of all-to-all 
propagators are important ingredients of our approach. 
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1 Introduction 



In recent years, there has been significant progress in determining the spectrum of 
hadrons containing a b quark, both experimentally [l}|5], and on the lattice [6]-|10|. A 
comparison of theory and experiment is of considerable interest for these hydrogen like 
systems, in particular since Heavy Quark Effective Theory (HQET) [Il| - [14| is applicable 
and is at the same time an important theoretical tool to isolate new physics in the flavor 
sector 
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The extraction of information on excited states from lattice simulations is a difficult 
problem, since the excited states appear as subleading exponentially decaying contri- 
butions to the lattice correlators, which at large times are dominated by the ground 
state and at small times by the combined contributions of arbitrarily highly excited 



states, and on top of this are affected by noise. A number of different methods 16 - 19 



have been proposed for overcoming this challenge; recently we have shown that with an 
efficient use of the generalized eigenvalue problem (GEVP), rigorous statements about 



the systematic error caused by the admixture of other states can be proven 20 . In 
particular, we have shown that for a suitable choice of Euclidean times t and to, the 
systematic error decays exponentially with an exponent given by an energy gap that 
can be made large by an appropriately chosen variational basis. 

Since their Compton wavelength is much shorter than any realistically achievable 
lattice spacing, b quarks cannot be simulated as relativistic quarks on a lattice. A de- 
scription of heavy-light mesons that is suitable for use in the context of lattice QCD sim- 
ulations is given by HQET with non-perturbatively determined parameters 21 22 . The 
parameters necessary to match HQET to QCD at order 0(l/mb) have been determined 
in the quenched approximation by our collaboration [23], and the non-perturbative de- 
termination for QCD with Nf = 2 dynamical quark flavors is far advanced [24| . 

In this paper, we present our results for the heavy-light meson spectrum from 
HQET using the GEVP method. In section [2| we give a brief review of our methods. 
Details of the simulations and data analysis procedures are given in section [3j and our 
results for the quenched heavy-light spectrum are presented. Section [4] contains our 
conclusions. Some technical details of our implementation of all-to-all propagators are 
relegated to the appendix. 



2 Methodological Background 
2.1 Non-perturbative HQET 

HQET on the lattice offers a theoretically rigorous approach to the physics of B-mesons 
since it is based on a strict expansion of QCD correlation functions in powers of 1/m-b 
around the limit mb — ?■ oo. Subleading effects are described by insertions of higher 
dimensional operators whose coupling constants are formally 0(l/mb) to the appro- 
priate power. This means that HQET can be renormalized and matched to QCD in a 
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completely non-perturbative way (25] , implying the existence of the continuum limit at 
any fixed order in the l/my, expansion. 

To fix the notation we write the HQET action at 0(l/mb) as 

^HQET = a^^{-2'stat(a;) - WkinOkin(a;) -WspmOspin(a:)} , (2.1) 

X 

^st.t{x) = ^iXx){Do + 6m)tPk{x), (2.2) 
Ckin(2;) = Vh(a;)D^^h(2;) , C'spin(x) = 'ilJh{x)a • Biph{x) , (2.3) 

where V'h satisfies ^-^tph = V'h- The parameters Wkin and cogpin are of order 1/mb, 
and 6m is the counter-term absorbing the power-divergences of the static quark self 
energy. The signal-to-noise ratio of large-distance correlation functions is significantly 
improved by replacing the link U{x,0) in the backward covariant derivative Dof{x) = 
[f{x) — U"^{x — aO,0)/(x — aO)]/a, with a smeared link [26 . 



Exponentiating the action of eq. (2.1) would give (non-renormalizable) NRQCD 



|27| ; in order to retain the renormaliz ability of the static theory, we treat the theory in 
a strict expansion in 1/mb, where the 0(l/mb) parts of the action appear as insertions 
in correlation functions. For the expectation value of some operator O this means 
(ignoring the possibility of explicit 1/mb operator corrections, which do not affect the 
energy levels [23| ) 

(O) = (O) Stat + ^kin « ^(OOkm(2;))stat + ^spin 0" ^(OOspin(a;))stat (2.4) 

X X 

where (O)stat denotes the expectation value in the static approximation. 

To fully specify HQET, the parameters 6m, cukm and Wgpin must be determined by 
matching to QCD. In order to retain the asymptotic convergence in 1/mb this matching 
must be done non-perturbatively, since for a perturbative matching at loop order /, the 
0((7^') truncation error of the static term is much larger than the power corrections to 
the static limit: 

epert(0«5''(mb)oc [26olog(mb/AQCD)]"' > • (2.5) 

mb-s>oo mb 

A fully non-perturbative determination of the parameters of HQET has been car- 
ried out in j23) . Here we employ the same discretization of QCD and HQET and in 
particular the determined values of Wkin and Wspin- For further details of the matching 



and discretization, the reader is referred to 23 



2.2 The Generalized Eigenvalue Problem 

In this section, we recall the relevant contents of [20]. Starting from some fields Oi{x) 
localised on a time slice and their momentum zero projection Oi{x) = Oj(xo), a 
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matrix of Euclidean space correlation functions, 

oo 

a,{t) = (O,(t)O*(0)) = ^e-^"*V'™Ci, iJ = l,...,N (2.6) 

n=l 

'Ipni = {lpn)i = {O\0i\n) En < En+l , 

provides the basis for the GEVP 

C{t)Vn{t,to) = Xn{t,to)C{to)Vn{t,to), n = 1, . . . , N , t > Iq . (2.7) 

Effective energies for the n-th energy level are given by 

E fit, to) = -dtlogXn{t,to) = -- [logXn{t + a,to)-logXn{t,to)], (2.8) 

a 

where a is the lattice spacing. Provided that to > the effective energies converge 
to the exact energy levels as |20l 



Ef{t, to) = En + 0(e-^^^+i." *) , AEm,n = - En. (2.9) 
In HQET, all correlation functions 

Cij{t) = C^^*(t)+L.kinC'^"(i)+^spinCf'^(t) + 0(^2) (2.10) 

are computed in an expansion in a small parameter, uj oc 1/m-b. Correspondingly, the 
energy levels expand as 

Ef{t,to) = i?f'^*^*(t,to)+a;xi^f'^(t,to) + 0(a;2) (2.11) 
ii;«(Mo) = a-' log .SuT'l^ (2-12) 

TpeS,^(+ + \ _ K.jt-.'to) A^(t + a, to) 

- Ar(Mo) " Ar(t + a,to) ^'-''^ 

where x G {kin, spin}, with the behavior at large time t < 2to, 

Ef^^^^\t, to) = + e-^^Nti,n * + ..., (2.14) 

Ef^-{t,to) = EI + [PI - /3r*t Ai?j^+i,Je-^^-+V"* + .... (2.15) 

where 

C''^\t)vT\t,to) = XT\t,to)C''^\to)vt''\t,to), (2.16) 



Kiji^'to) 
Xii^\t,to) 



vT\tM) , [[>^t\tMr c-{t) - c-{to)]vT\t,to)) . 



with {v^^''\t, to) , C"'*^*(to) v^^'^^t, to)) = Smn- We note that the GEVP is only ever solved 
for the static correlator matrices. 
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2.3 Mass splittings to 0(l/mb) 



With the static Lagrangian eq. (2.2), all HQET energies satisfy exactly En = En\g,^^Q + 
^ log(l + aSm), and the power divergent dm drops out in energy differences. Since we 
only consider these in this paper, we need just and cjgpin of [23) . 

The excitation energies at the static order of HQET are given simply by the differ- 
ences of the static energies, 

A£;^*f = Ef^^^ - Ef"^ , (2.17) 
and at order l/m-b, the excitation energies of pseudoscalar Bs states become 

^^HQET ^ ^^stat _ ^stat) ^ ^^.^^E^in _ ^kin) ^ ^^^^^(i^spin _ ^spin^ _ ^2.18) 

At the static order, the masses of pseudoscalar and vector states are degenerate due 



to spin symmetry 13 . This degeneracy is lifted at the 0(l/m]-,) level by the contribution 



from Ospiru giving & Bg — B* mass difference of 

3"^spm-'^l 



A£;p_v = ta;spm^r° • (2-19) 



3 Simulation details and results 
3.1 Lattice parameters 

We use three quenched ensembles of 100 configurations each, generated using the Wilson 
gauge action with parameters /3 = 6.0219, 6.2885 and 6.4956. The physical volume was 
kept constant at L ~ 1.5 fm, leading to L/a = 16, 24, 32 for the three ensembles. We 
used time extent T = 2L throughout. 

The static quark is discretized on each ensemble with both the HYPl and HYP2 



[26,28,29 actions. The light valence quark is discretized by a non-perturbatively 0(a)- 
improved Wilson action [30|31| , and its mass was fixed to the strange quark mass, giving 
Ks = 0.133849, 0.1349798, 0.1350299, respectively [32]. A summary of the simulation 



parameters used, with the corresponding values of Wkin and Wspim is given in Table [T| 



3.2 Measurements of correlation functions 

The strange quark propagators are computed through a variant of the Dublin all-to-all 
method |34j. We use approximate instead of exact low modes (the method remains 
exact) and employ even-odd preconditioning in order to reduce the size of the stochas- 
tically estimated inverse of the Dirac operator by a factor of 2. The reader is referred 
to appendix [K\ for details. 

The interpolating fields are constructed using quark bilinears 

Ok{x) = ^^{x)jol5i^l'\x) (3.1) 
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/3 


ro/a 


{L/af X T/a 


Kg 


Nl 






Wspin / O 


6.0219 


5.6 


16^ X 32 


0.133849 


50 


2 


0.330(5) 


0.55(4) 


6.2885 


8.4 


24^ X 48 


0.1349798 


50 


2 


0.439(7) 


0.71(5) 


6.4956 


11.0 


32^ X 64 


0.1350299 





4 


0.553(9) 


0.87(6) 



Table 1: Parameters of the simulations: inverse coupling /3, approximate value of the 
scale parameter tq |33] in lattice units, spacetime volume, hopping parameter for the 
strange quark mass 32 , number of low-lying eigenmodes and number of noise sources 



used in the all-to-all 34 estimate of the strange quark propagators, and values of the 
HQET couplings Wkm and Wspm 23 for the HYP2 action in lattice units. 



built from the static quark field ^i^hix) and different levels of Gaussian smearing [35] for 
the light quark field 

V'f'^(x) = {l + KGa^ ^f'Ux): (3.2) 
where the gauge fields in the covariant Laplacian A are first smeared with 3 iterations 



of (spatial) APE smearing [36,37 to reduce noise. At /3 = 6.2885, we use Rk = 0, 22, 
45, 67, 90, 135, 180, 225 with kq = 0.1. At the other values of /3, we rescale the values 
of Rk used so that the physical size ?^phys,fc ~ 2a\/KG-Rfc of the wavefunctions is kept 
fixed; in particular we keep rmax = ^phys,7 ~ 0.6 fm. 

For these bilinears, we compute the following correlators: 



Cf'it) = Y.{0^{xo + t,y)0]{x))^^^^, 

(3.3) 

.Wspm^^^ = (O^(^0 + t,y)O-(x)Okin/spin(^)>3tat 



with the 0(l/mb) fields defined in eq. (2.3). 
3.3 Determination of energies 



The correlator matrices of eq. (3.3) are "thinned" to form a sequence of N x N matrices, 
where N G {2,..., 7}, by selecting only entries from a certain subset X^v of indices 
where Xjv = {1,7}, {1,4,7}, {1,3,5,7}, {1,2,4,6,7}, {1,2,3,5,6,7}, {1,2,3,4,5,6,7}. 
An alternative procedure, used in [20] is "pruning" [38], where the N x N matrices are 
formed by projection onto the subspaces spanned by the lowest eigenvectors of C{ti) 
for some small tj. We decided not to use this version since it introduces a dependence 
on the relative normalization of the fields Oj. Moreover, with thinning we found a 
somewhat faster convergence to the plateau for the ground state energy compared to 
the pruning version [39] when the normalization of Oi in l20] is used. We note that the 



fact that we found the same low-lying energy levels with thinning as with pruning is a 
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n 


2 3 


4 5 6 


roK^n - ) 


1.50(5) 2.7(1) 


4.0 5.0 6.0 



Table 2: Rough estimate of the static spectrum. For n < 3, the gaps come from the 
continuum hmit; for n > 4, the rough estimate used to stabihze the fit is quoted. 

further confirmation that the GEVP is quite robust against changes of the variational 
basis employed. 

For each of the resulting N x N correlator matrices, we solve the static GEVP and 
compute the static and 0(l/mb) energies. This gives a series of estimates E^'^^'^^{N, t, to), 
E^'^^'^{N,t,to) and E^'^^^^{N,t,to) with associated statistical errors, which we deter- 
mine by a full Jackknife analysis. 

To arrive at final numbers for En we also need to estimate the size of the systematic 
errors coming from the higher excited states. To do this, we first perform a fit of the 
form 

to the GEVP results for Ef'''^''\N, t, to), fitting the data at iV = 3, . . . , 5, ^ro < to < 6a, 
^0 ^ ^ ^ 2to and n = 1, . . . , 6 simultaneously. The stability of the fit is enhanced in the 
following manner: First we perform an unconstrained fit to extract E^^^ —Ef'^^ for n = 
4, 5, 6 for each lattice spacing and action and compute a rough average of ro^E^^^—Ef'^^) 
for these values of n. Then we repeat the fit, constraining tq^E^^^ — Ef"^^) for n > 4 
to the previous average (this renders the fit linear). For n < 3 this is unnecessary as 
these levels are well determined at each individual lattice spacing. We list the extracted 
values of ro(£;^*^* - Sf^*) in Table 

Finally, using the values of E"^*^* and determined from this fit as (fixed) input 

parameters, we fit E^'^^'^{N,t,to) and E^'^^^^{N,t,to) by 

Ef^^'-{N,t,to) = El'- + e:i^^^it) (3.5) 



E^''' + 



okin ostat j. / i^kin irikinN 



( iT'stat 77stat\ + 



Kf'^P'"(iV,t,to) = ^r"" + en (3-6) 



ET'" + 



ospin ostat j. / irispin pspm\ 
Pn,N ~ Pn,N \-^N+l ~ ) 



/ TT^stat rpstat\4 



in the same manner. 

While the fitted results are quite stable, we consider them as rough estimates only, 
since our fits include only the leading exponential correction, and there are systematic 
effects from higher corrections to the GEVP. We therefore employ the fits only to esti- 
mate the size of the leading corrections. For a reliable estimate of the energy levels, we 
calculate plateau averages from t = tmin ^ on at each N and to, and take as our final 
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Figure 1: Illustration of some plateaux. Left: £1^^*^* (bottom half) and E^^^^^ (upper 
half). Here N = 5 and within each group the lattice spacing is decreasing from top 
to bottom. Dotted lines represent the global fit, while dashed lines indicate the cho- 
sen plateau. Right: the ground state spin splitting. The ordinates of the points are 
shifted by an arbitrary /3-dependent offset in order to make the different plateaux visible 
separately. 



estimate that plateau average for which the sum utot = Cstat + Csys of the statistical error 
Cstat of the plateau average and the maximum systematic error dsys = e(imin) becomes 
minimal, subject to the constraint that Usys < gCstat- We impose the latter constraint 
in order to ensure that the total error is dominated by statistical errors. 

An illustration of the more problematic plateaux is shown in Fig. [T| It is rather 
clear that without some analysis of corrections due to excited states it is very difficult 
to locate a safe plateau region at least for n = 3. 

Our results are given in Table[3) besides the plateaux found by the method described 
in the preceding paragraph, we also show the results of the global fit, which in almost 
all cases agrees very well with the final plateau value. 

3.4 Continuum limit 

We now turn to the continuum extrapolation of the level splittings. Using the fact 
that the static actions employed are discretizations of the same continuum theory, we 
perform a combined continuum limit by fitting a function of the form (k = 1,2 for 
HYPl, HYP2 actions) 

^i,fc(aAo) = + Q,fc(a/ro)'' (3.7) 
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HYPl HYP2 



p 


WUot;! VdUiC 


Fit 




Fit 

-Tit 


JT IdLCdU. 




^ IT'S tat 














0.708(2) 

0.893(5) 


0.711(5) 

0.92(2) 


0.678(2) 

0.868(5) 


0.680(6) 

0.89(2) 




a^Ef"" 

„2 pkin 
2 


0.743(1) 


0.740(2) 


0.774(1) 


0.771(2) 

U.OOl ± I 




2 pspin 


-U.UOUUl (J J 


-yj .yjziijoy t j 


-yj .yj^ooy^ j 






2 iTispin 
a ^2 










6.2885 


aEf^^ 
aEf^^ 


0.3319(2) 

0.507(1) 
0.648(3) 


0.3328(7) 

0.515(4) 
0.67(1) 


0.3032(2) 

0.478(1) 
0.614(3) 


0.3041(7) 

0.485(4) 
0.625(7) 






0.6479(5) 
0.682(4) 


0.6481(4) 
0.68(1) 


0.6743(4) 
0.704(4) 


0.6745(3) 
0.70(1) 






-0.0127(1) 


-0.0126(2) 


-0.0129(1) 


-0.0130(3) 






-0.0115(9) 


-0.011(1) 


-0.011(1) 


-0.012(1) 


6.4956 




0.2742(3) 


0.275(1) 


0.2482(3) 


0.249(1) 




^^stat 

a^Ef"" 


0.409(2) 
0.518(3) 
0.5999(5) 

0.620(3) 


0.405(8) 

0.52(2) 

0.5997(3) 

0.625(6) 


0.384(2) 
0.491(3) 
0.6240(3) 

0.645(2) 


0.381(8) 

0.50(2) 

0.6229(7) 

0.645(4) 




a'^Ef' 


-0.0081(1) 


-0.0079(5) 


-0.0076(1) 


-0.0074(4) 




a^Ef" 


-0.0108(9) 


-0.005(4) 


-0.0098(8) 


-0.007(2) 



Table 3: The measured values of the HQET energies in lattice units. Shown are both the 
values obtained from a global fit ("Fit") and from our more conservative plateau selec- 
tion ("Plateau") as described in the text, for both the HYPl and HYP2 discretization 
of the heavy quark. 
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/3 


= 6.0219 


/3 = 6.2885 


/3 = 6.4956 


cont. limit 




HYPl 
HYP2 


1.50(3) 
1.51(3) 


1.52(9) 
1.51(3) 


1.46(8) 
1.47(8) 


1.50(5) 


roAEff 


HYPl 
HYP2 


2.67(10) 
2.72(10) 


2.78(9) 
2.68(6) 


2.7(2) 
2.8(2) 


2.7(1) 




HYPl 
HYP2 


0.20(2) 
0.21(2) 


0.14(4) 
0.08(4) 


0.17(5) 
0.13(3) 


0.03(6) 


roAE^^^^ 


HYPl 
HYP2 


1.70(4) 
1.72(4) 


1.66(6) 
1.59(6) 


1.63(10) 
1.59(9) 


1.56(8) 


^oAi?2jl 1 continuum ~^ 


continuum 


1.54(9) 


roA£'p_v 


HYPl 
HYP2 


-0.093(7) 
-0.114(9) 


-0.083(6) 
-0.103(8) 


-0.089(8) 
-0.092(7) 


-0.075(8) 



Table 4: The energy level differences in units of the scale rg. Shown are the results at 
each /3 for both static-quark actions, together with their common continuum limit. 

to our dimensionless quantitites $i G {roAE^I^f , roA£'|*f , roA£;p_v, roAS^^^'^}. Since 
0(a)-improvement is fully implemented in the static approximation, we use powers 
si = S2 = 2. On the other hand, the l/rrib corrections have 0(a) discretization errors, 
yielding S3 = 1 for the observable A£^p_v. For A^^^^^""", there are two possible ways of 
taking the continuum limit. First we extrapolate AE^f" and the 0(l/mb) contribution 
separately to a — t- and add the continuum limits afterwards. Here we set sl*^* = 2 
for the static part and = 1 for the 0(l/mb) correction. Second, one can form the 
combined AE^^^ at each (3 and take the continuum limit of the combination. The 
linear term in a which is present in the combined data, is suppressed by a factor 1 / mb ■ 
Given in addition the flatness of the data in a (see Fig. [2]) we just use S4 = 2 for the 
combination, assuming that this term dominates. 

Figure [2] shows the approach to the continuum limit for the static and full HQET 
energy splittings. The results for HYPl and HYP2 are distinguished by color; the static 
splittings are shown as circles, whereas the full HQET splitting AE2 i is shown as 
diamonds. Also shown are the fits to the continuum limit together with their error 
bands. We see that the approach to the continuum limit is rather flat in particular for 
A£^2,i) and that the 0(l/mb) corrections constitute only a small shift of the energy 
splitting between the first excited and ground states. In Fig. [3] we show the approach 
to the continuum limit for the spin splitting A£'p_v in the same fashion. 

3.5 Results 

Our findings for the continuum values of the level splittings are summarized in Table |4] 
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k i 
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-f=— i — 


a 


— HYP1 

— HYP2 



'O.OOO 0.002 0.004 0.006 0.008 
a2 [fm2] 



Figure 2: Plot of the continuum limits (stars) of AEf^'f^ (circles) and i (dia- 
monds). Shown are the results for both HYPl (red, shifted to the left) and HYP2 
(blue, shifted to the right). 



For the static spectrum, we obtain 

roAEff = 1.50(5) (3.8) 
roAEff = 2.7(1), (3.9) 

corresponding (using ro = 0.5 fm) to AEff = 594(21) MeV and AE^^f = 1076(48) 
MeV in good agreement with the results of |10| at a fixed lattice spacing. These numbers 
in physical units are meant as a rough illustration, since no quenching error is attached 
to them. 

At 0(l/mb), we obtain tqAE^^^"^ = 1.54(9) when taking separate continuum 
limits for the static and 0(l/?7ib) energy differences, and 

roAEf^^^ = 1.56(8) (3.10) 

when combining static and 0(l/mb) energy levels before taking the continuum limit. 
The results from both procedures agree within errors, indicating that the 0(o) term 
omitted in the combined extrapolation is at most a minor source of systematic error. In 
physical units (again using ro = 0.5 fm), our results correspond to AE^^^'^ = 606(35) 
MeV and 617(31) MeV, respectively. 

For the Bg — B* mass difference, we find 

roAEp_y = -0.075(8) (3.11) 
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Figure 3: Plot of the continuum limit (star) of A£'p_v for HYPl (red, shifted to the 
left) and HYP2 (blue, shifted to the right). 

giving A£^p_v = —29.8(3.2) MeV via tq = 0.5 fm. This is to be contrasted to the 
experimental value 40 of —it^b*= —49.0(1.5) MeV. We note that while 0(l/m^) 



effects may not be entirely negligible for such a small splitting, the difference is too large 
to be entirely attributed to these. Instead, a genuine quenching effect is involved. 

The spin splitting between the first radial excitations in the pseudoscalar and vector 
channels is 

roAEp>_y> = -0.056(27) , (3.12) 

which is compatible with the spin splitting between the ground states. 

Our results confirm the expectation that the 0(l/mb) contributions are small 
compared to the static results. To get an idea of the size of higher order correc- 
tions in 1/mb, we have considered the dependence of the energy level splittings on 



the matching conditions chosen in the determination 23 of the HQET parameters cj^in 
and Wspin- We find that the dependence is very minor, with a maximum deviation of 
6\rQAEp-Y\ = 0.005, less than the statistical errors, and 5\roAEl^^\ = 0.002, much 
less than the statistical errors. Comparing this to the naive power counting estimate of 
|0(l/m^)| ~ (l/(romb)^ ~ 1/100, we see that the tested 1/m^ terms are as small as ex- 
pected or even smaller. Note that 1 /ro ~ 400 MeV is indeed a typical non-perturbative 
QCD scale. 



11 



4 Conclusions 



In this paper, we have reported results from quenched lattice QCD for the spectroscopy 
of the low- lying excited states of the Bs and B* systems. An application of the gen- 
eralized eigenvalue method with all-to-all propagators to non-perturbative HQET at 
0(l/mb) allows us to extract precise results for the energies of the lowest-lying radial 
excitations as well as for the Bg — B* splitting. However, we emphasize again that a 
careful analysis of systematic errors due to excited state contaminations is necessary. 

A first relevant observation to be pointed out concerns the renormalizability of 
HQET. Unlike for QCD on and off the lattice, there is no proof of renormalizability 
of the theory to all orders of perturbation theory. However, we find that in our non- 
perturbative computations the divergences cancel after proper renormalization of HQET 
p3] . The left over lattice-spacing dependence in Fig.[2| Fig. |3]is very flat. To appreciate 
this, note that 

r^E^^ w (24 , 47 , 77) for a = (0.1 , 0.08 , 0.05) fm (4.1) 

as seen in Table [s] and weaker but still very prominent divergences are present in E^'^^ . 
In other words we find strong numerical evidence for the renormalizability of the theory; 
in fact also the universality of the continuum limit is demonstrated in the figures. It is 
also worth emphasizing that the present demonstration is the first time the continuum 
limit is taken for mass splittings in HQET. 

We find the physical 0(l/mb) corrections to be small throughout. 

The precision attained, in particular when taken together with the relative small- 
ness of the 0(l/mb) effects, indicates that non-perturbative HQET combined with the 
use of the GEVP for data analysis is a reliable method for determining B meson spectra. 
We intend to apply it to the Nf = 2 case in the near future. In this context one should 
remark that we were able to achieve good precision using only 100 configurations in our 
quenched study. Therefore we do expect to be able to decrease the errors for dynamical 



fermions. However the influence of topological modes being updated only slowly 41 
needs to be controlled or better algorithms with a faster decorrelation need to be used. 
A promising proposal has been made in i42j. 
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A All-to-all propagators 



In this appendix, we explain the details of our implementation of all-to-all propagators, 
which follows the idea of 134] with some useful improvements. 



A.l Even-odd preconditioning 

To reduce the computational effort and storage requirement for the matrix inversions, 
we consider even-odd preconditioning of the (hermitian) Wilson-Dirac operator Q = 
2k^5D. With even/odd ordering of the sites one has a block structure 

( Mee Meo 
^ V Moe Moo 

where Mee (Moo) differs from unity by the clover term on the even (odd) sites, and Mq 
(Meo) is the hopping term. Defining 



B 



le -M-^Meo 



lo 

the preconditioned matrix B'^QB is block-diagonal and the propagator can be factorized 



as 

where Qee 

= 75Mee is diagonal in space-time, and Qoo = lb{Moo- MoeM^J" M^o) = QIo- 

A.2 Approximate low modes and a stochastic estimator 

We consider an orthonormal basis {vi : i = 1 . . . N^} of an A'^^ dimensional subspace 
{Nl > 0) of all fermion fields which live only on odd sites. Defining the projectors 

Nl 

^i^ = ^Vi-v\ and Ph = 1o-Pl, 

i=l 

we can write 

Nl 

Qoo = Qoo'(Pl + Ph) = Y.^Qo}h) ■ v\ + Q^o'Ph . (A.2) 

i=l 

A natural choice for ■Oj are approximate eigenvectors of the low-lying eigenvalues of Qoo 

QooVi = hvi + r-i . (A.3) 
with \\vi\\ = 1 and ujrfc = 0. Then, the part Qoo^Pl in (A.2) is expected to approximate 



the long-distance behaviour of the propagator 34 ,43 , and the inversion^ needed in 
(QooVi) are cheap. 



^ Since we do not explicitly use Qoaiii ~ ^Vi, the errors ||fi|| in ( A.3I are allowed to be large. In 



practice, we require ||ri|| < 0.001 • |Aij, and take only as start vectors for the inversion. 
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On the other hand, we can introduce a stochastic estimator for Ph- We take 
random vectors r]i with 



{'ni,a ffj,p)v 



0, 

0, 



(A.4) 
(A.5) 
(A.6) 



where a, (3 denote combined (color, Dirac, and site) indices, and (.)^ is the average over 
ry. The relations (A.4)-(A.6) hold, for instance, in the case of a U(l) noise 

?7i,Q, = exp(i(/)i,Q,) , 

where (jji^a are independently uniformly distributed in [0, 27r) (while for Z2 noise the 
average in (A.6) would not vanish for i = j). Thus, the second term in (A. 2) can be 
written as 



Qoo Ph = ^ 5] {Qoo' Ph m ■ vl)r 



(A.7) 



i=l 



and the estimator of Q^^ can be written as a sum of dyadic products 



Qoo= E 



(A.8) 



witl 



Wi 

m 



Qoo ! 

Qoo Ph Ui , 



Ui 
Ui 



Vi 



{i = l,...,NL) 



The full propagator is then obtained from (A.l). Since B connects only adjacent 
time slices, the block does not contribute to the propagator between sites with time 
separation |xo — yo| > 2a. In this case, we can simply write 



Q 



-1 



E < 

1=1 



Wi ■ u. 



(A. 



with 



Wi 



B 





Wi 



and Ui 



B 





Ui 



Even-odd preconditioning can be seen as a form of dilution 34 , since there are 
only half as many components of the noise field r] in the even-odd preconditioned case 
as without preconditioning. Note, however, that unlike other dilution schemes, even-odd 
preconditioning does not increase the number of inversions needed. 



^ One may also use Wi 
not tested this option. 



Qoo and iii 



PiiVi for i ^ Nl + 1, 



, Nl + N„ , but we have 
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A. 3 Time dilution 



In addition, we use the more conventional time dilution scheme. It is implemented by 



replacing Ph in (A. 2 ) by Ph Ylt -^t where Pf projects on the components corresponding 
to (odd) sites with time coordinate t. Then, an independent stochastic estimator is 
introduced for each term 



V 



i=l 

where the noise vectors r]ti have non-vanishing components only for (odd) sites on time- 
slice t. 



Note that due to the hopping term in B the full propagator (A. 9) from time slice 
xq to uq receives contributions from noise vectors r]ti on three time slices, t = xq, xq it a, 
i.e. three inversions are required for the propagator from one time slice xq. However, a 
total of T inversions is sufficient, and hence no extra effort is required, if one computes 
the propagator for all xq, as we do in our measurements. 

Analysing the variance of a heavy-light two-point correlator as described in [44], 
one sees that the variance with time dilution decays roug hly as e-^^o-^^o)™-, while the 
expression without time dilution contains pieces independent of xq — i/q. This renders 
time dilution very profitable. 
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